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Abstract 

Slater type orbitals were used to construct the overlap and the Hamiltonian core matrices; we also found the values of 
the bi-electron repulsion integrals. The Hartree Fock Roothaan approximation process starts with setting an initial 
guess value for the elements of the density matrix; with these matrices we constructed the initial Fock matrix. The 
Mathematica software was used to program the matrix diagonalization process from the overlap and Hamiltonian 
core matrices and to make the recursion loop of the density matrix. Finally we obtained a value for the ground state 
energy of helium. 
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Introduction 

The Hartree Fock method (Hartree, 1957) also known as Self Consistent Field (SCF) could be 
made with two types of spin orbital functions, the Slater Type Orbitals (STO) and the Gaussian 
Type Orbital (GTO), the election is made in function of the integrals of spin orbitals implicated. 
We used the STO to construct the wave function associated to the helium atom electrons. The 
application of the quantum mechanical theory was performed using a numerical method known 
as the finite element method, this gives the possibility to use the computer for programming the 
calculations quickly and efficiently (Becker, 1988; Brenner, 2008; Thomee, 2003; Reddy, 2004). 

It is known that the effects of electron spin produces an antisymmetric wave function and 
that the charge of the nuclei is screened by the other electrons, then the effective charge that 
affect the interactions electron-nuclei and electron-electron is all of this effects are considered 
on the Hartree Fock Roothaan approximation. Dealing with Restricted Closed shell 
approximations means that we are considering only an even number of electrons, with all 
electrons paired such that n=N/2 spatial orbitals are doubly occupied (Szabo, 1996; Levine, 2008; 
Bransden, 2003). 

Hartree Fock Roothaan equations 

The initial point of view of this approach is the independent particle model, according to which 
each electron moves in an effective potential which takes into account the attraction of the 
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nucleus and the average effect of the repulsive interaction due to the others electrons and their 
spin (Hartree, 1957; Rosetti, 1974; Szabo, 1996; Levine, 2008; Bransden, 2003). 

The general procedure for finding a wave function of order zero (ground state) of a set of 
electrons, is to multiply the hydrogen-like wave function of each one electron 

— f n{ r nA>9n) 

Where the hydrogen-like orbitals are 




Where R nl ( r ) is the radial function and Y" 1 (6,(p) are spherical harmonics. This wave function 
has total nuclear charge Ze, with Z number of protons and e the charge of an electron. 

Each electron in a many electron system is described by its own wave function (Bransden, 
2003). According with this approach the trial wave function of the A-electrons is a Slater 
determinant ®, which is a function of spatial orbitals and spin states. The central modification 
made at this point is to assume that an electron in a multi-electron atom not perceive the whole 
nuclear charge, as this is partially shielded by the other electrons. Therefore it is proposed an 
antisymmetric wave function of the form (Szabo, 1996; Levine, 2008; Bransden, 2003) 

Where g n represent the spin orbital, but now depends of shielded nucleus charge £, which can be 
expressed as a Linear combination of Atomic Orbitals (LCAO), depending on the base functions 
( x ) as follows 


(!) 

/=i 

Where the base function / are a complete set of Slater Type Orbitals (STO), which is defined in 
atomic units as 


Z = 



( 2 ) 


in this STO n represents the principal quantum number of the electron orbital for which the 
function / is defined, and C, is an effective electron charge (Rosetti, 1974). 

The Hartree-Fock-Roothaan equations based on the STO functions are 


b b 


f =H are +yyp 

rs rs / j / ^ tu 


t =1 i/=l 


r^tu)- 


( 3 ) 
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in equation (3) F rs represents a general element of the Fock Matrix, (rs\tu) are the electronic 
repulsion integrals and P, u are the elements of the density matrix, that are written as 

nil 

P,u=^L c Vi ; t = ; u = 1,2,. ..,6 (4) 

./= 

electronic repulsion due to the charge and spin of the electrons expressed in the spin orbital x are 
given by the integrals 




JJ Zr( 1 )z.( l )A 2 )x»( 2 ) dZ ' dT 


r u 


( 5 ) 


These integrals represent the expected value of the interaction l/r 12 relative to the state 
X r (i)X s (j) such that the electron i is in spin-orbital x r and the electron j is in spin-orbital % s . The 
H""' e are the elements of the Hamiltonian operator that are given by 


ht=( x,y)\H core \xy)) 


( 6 ) 


Where the Hamiltonian operator written in function of the screened charge £ is 


H OTO 



(V) 


The Roothaan eigenvalue equation is then expressed as 

T c A F rs - V‘ S 'J = ° => Determinant(F rs -s i S rs ) = 0 (8) 

s=l 


Finally the elements of the overlapping matrix are the integrals given by 

S rs= {z r \X s )=\z* r {rA<p)z s {rA<p)dT 


( 9 ) 


Integral evaluation 

There are two kinds of integrals, the one-electron as are the overlapping integral S rs and the 
Hamiltonian core H c r ° re , and the two-electron given by the electron repulsion integral (rs\tu) , each 
one has different treatment. 

For the one-electron integrals, we used the formula shown below, that can be obtained 
integrating by parts 


J x h e- qx dx 
0 


b\ 


( 10 ) 
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and for two-electron integral we used (10) and additionally these two ones 


x 2 e~ bx dx = 


r x 2 2x 2 A 

- 1 - TT H - T 

yb b 2 b 3 J 


( 11 ) 



( 12 ) 


Since in the application of equation (5) to the spin orbital of two electrons recursively we 
have the integral that is shown below 



( 13 ) 


This integral can be transformed to an evaluation formula as is described next. 

We know that r n represents the distance between the two electrons and that dr 2 is the 
differential radial in the space of the electron 2, then we have 


A = \e~ ar 'r 2 dr x B (14) 

0 


Where B is given by 


B- 


1 ri 

-jr 2 v 

r J 


2 dr 2 + 1 1 


2 dr 


Then applying (\ref {11}) to the first part of B and (12) to the second part, we have 



b 2 b 3 r { 


(15) 


Substituting (15) in (14) and simplifying we obtain 


A 


771 e ~ ar ' r \d r \ --jr]e~ (a+b)r 'r 2 dr l - e~ (a+b)r 'r,dr x 

O 0 D 0 0 0 


(16) 


Using equation (10) to evaluate (16), finally obtain a formula for the electron repulsion integral, 
as 


_2_2_ 2 

b 3 a 2 (a + b) 3 b 2 (a + tif b 3 


(17) 


4 














EJPE 


European J of Physics Education Volume 5 Issue 1 2014 Acosta, Tapia, Cab 


Therefore, if we want to calculate the ground state of helium atom with the SCF procedure, 
using the orbitals exponent £ =1.45 y £,=2.91 (Roetti, 1974), the spin orbital functions are, 
according with equation (2) 


Xl =2S}e<'X 


Z2 =2£V^F 0 ° 


(18) 

(19) 


Applying (5) to the state (l l|l l), we have 


„i„ )= jf4MMbki(9 df|dfi 

^12 


Substituting spin orbital and simplifying 


(ll|ll) = 16£ 6 jV 24ri,i 



( 20 ) 


with dr ] = f SinOfrfdf^ , dr 2 = r 2 Sin0 2 dr 2 d0 2 d(f> 2 and 


n 2 2;r 

J(7 0 °) sin OdOj d</> = \ 
0 0 


Equation (20) is written as 


(ll|ll) = 16£ 6 J e-K'rfdrJ 
0 0 




r 2 dr 2 


Then using (17) with a = b = 2£ 


(ll|ll) = |£=0.9062 

for the state h.2\22 \, we have the same result, but with £ 

(22|22) = |<£ =1.8188 

Again applying (5) to the state (l2|l2), we have 


( 21 ) 


( 22 ) 
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Substituting spin orbitals and simplifying 


°° °° /=> _ (^l +l ^2) r 2 

(12|12) = 16 f l 3 £je- (c ' +C2h r l 2 dr l j - r 2 dr 2 


with a = b = (g l +f 2 ) substituting in (17) and simplifying 


(12|12) = 2 °^ ^ = 0.9536 


(23) 


Once again applying (5) to the state (l \\ 2 l ), we have 


ll | 22 ) = jj 


Substituting spin orbitals z,, Xi and simplifying 


r 2 ^2*2 


11122) = \ey x y 2 J e- 2( ' r 'r 2 dr,\^—r 2 dr 2 


with a = 2£ x and b = l^ 2 substituting in (17) and simplifying 

1 1 1 


(ll|22) = £% 3 




= 1.1826 


(24) 


Again applying (5) to the state (l 1|12^, we have 


(i !| 12 )=j j iMMM 2 !*,*, 


Substituting spin orbitals Z\ > X 7 and simplifying 


9 300 co 

(l 1|12) = 16^ C 2 2 j e~ 2(,r ' r 2 dt\ J- r 2 dr 2 


with a = 2£i and b = g l +£ 2 substituting in (17) and simplifying 
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W\U) = \6QQ 


(3£ + £ ) 3 -4 (£ + QC[ - 4(3^ + QC[ 

2(3£i + 4*2 ) 3 (4”i + C 2 ) 3 Ci 


it's easy to observe that 


(l 1|12) = (l2|l l) = (l l|2l) = ( 2 1 |l l) = 0.9033 


Finally applying (5) to the state (l2|22), we have 

(I2| 22 )=jj 




( 25 ) 


(26) 


Substituting spin orbitals , % 2 and simplifying 


,“ 2 ^2 r 2 


(12122) = 16£ 2 £* J r y- {(l+(l)r ' dr, J-— r*dr 2 


with a = £ + £ 2 and b = 2^ 2 substituting in (17) and simplifying 


(12|22) = 16^| 


4 c 3 2 (c 1 + c 2 y 2 (^+ 3 c 2 yc 2 2 


it's easy to observe that 


(12|22) = (22|12) = (2l|22) = (22|2l) = 1.2980 


(27) 


The one-electron integrals are simpler than the two-electron integrals, then for the 
overlapping integrals, we have that S n =(z l |^) = l ; S 22 = {z 2 \x 2 ) = 1, and for S, 2 = {zi\x 2 ), utilizing 

(9) 

3 3 oo 

S n ={ Xl \ X2 ) = 4C}C}\e- (C ' +(& r 2 dr 


And evaluating with (10) it is obtain 


S l2 = 


(G+6,) 3 


= 0.8366 


( 28 ) 


The Hamiltonian core operator ( H core ), for helium atom, are obtained using (7) 
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w are =-iv 2 -- 

2 r 


Since each element of the Hamiltonian is calculated with (6) 


( 29 ) 


"r=U 


I v *_2 

2 r 




Where the Laplacian for this element is 


V 2 ( e -^) = - l T —\r 2 — («-*') 
V ’ r 1 dr dr y ' 


Then integral for the Hamiltonian core element H™ re is 




Simplifying and evaluating using (10), it is observed too that H™ re has the same expression 


HT = 2^1 - 2 ^,=--1-8488 


-2^ =-1.5860 


(30) 


Then the integral for the Hamiltonian core element H“ re is 


HT = 4&ll(e*) ^( 2 -Q)-- r (e^)rdr 


Which is evaluated using (10) 


Htr=4tf£. 


QC 2 - 2 £ - 2^ 2 

(£+^) 3 


H"[ e = -1.8826 = H™ 


(31) 


Matrix Method Description (MMD) 

As we can see the overlapping matrix (S) is not orthogonal, 

T 1 0.83661 

S = 

0.8366 1 


(32) 
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That is because the non-diagonal elements are different from cero. Therefore the procedure 
initiates with the diagonalization of matrix S, in this particular example the eigenvalues are 
S = {1.8366,0.1634} , here it is extremely important to note that the eigenvalues and their 

eigenvector associated, must be sorted from lowest to highest, this is to arrange the spin orbitals 
in the same sense. 


\ =0.1634 □ eigenvector ^ =( 0.707106781,0.707106781) 


Therefore, we build a unitary matrix (U) with the sorted eigenvectors taken as columns of U 


U = 


-0.707106781 0.707106781 
0.707106781 0.707106781 


What we have done is to find a basis, complete and orthogonal, but not orthonormal. If we 
want to facilitate all the operations involved in the SCF procedure, we must work with a unitary 
transformation matrix ( X ), then all the elements of matrix U must be divided by its corresponding 
eigenvalue’s squared root. 


Ji — 

U 2l U 22 

For the helium this matrix ( X) is 


- 1.749278568 0.521768326 
1.749278568 0.521768326 


(33) 


An additional matrix (D) is the one built with the eigenvalues in the principal diagonal 


D= 0.1634 0 

0 1.8366 

Some matricial proof we can make is that XDT = /, and UDU ' = S . 

All the values founded in previous section (Szabo, 1996, Levine, 2008, Bransden, 2003, 
Roetti, 1974), are part of the Fock matrix, then each element of this matrix are calculated using 

( 3 ) 


T7 — Tjcor 

Ai ~ M n 


+ ^(( 11 l 11 )-^( 11 l 11 )] + ^(( 11 l 12 )-^< 12 l 11 )] 

+ P 21 f(ll|2l)-I(ll|2l)l + J P 22 f(ll|22)-I(l2|2l) N 
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F u =- 1.8488 + 0.4531/], + 0.9033/] 2 + 0.7058P 22 (34) 

F n = H™ +^.f( 12 |ll) — \{ U |l2)] + /? 2 [(12112)-1(12|12)1 
+P 21 ^(l2|2l)-i(ll|22)j + P 22 ^(l2|22)-i(l2|22)j 

F n =- 1.8826 + 0.45165/],+0.8391/] 2 +0.6490P 22 = P 2 , (35) 

F 22 = H™ +^^(22|ll)-|(2l|l2)j + i>^(22|l2)-i(22|l2)j 

+P 21 ^(22|2l}-i(2l|22)j + P 22 ^(22|22)-I(22|22)j 

F 22 =- 1.5860 + 0.7058/],+1.2980/] 2 +0.9094P 22 (36) 

At this point we must obtain an initial guess at the density matrix P, which commonly are 
zeros or ones (we have taken ones), therefore the initial matrix F is 

[ 0.2134 0.057151 

F = 

0.05715 1.3272 


Transforming the matrix F to the orthonormal space using X, this is what we mean FI = X ] FX. 
Now we must diagonalize the Fock matrix FI, again eigenvalues and eigenvectors associated 
must be sorted from lowest to highest. 

s, = 0.202241 □ eigenvector v; = (-0.237268,-0.971444) 

s 2 - 4.61274 □ eigenvector ^ = (-0.971444, -0.237268) 

The matrix obtained named Cl, is in the orthonormal space 

[ 0.237268 -0.971444] 

Cl= (37) 

-0.971444 -0.237268 v 7 


Then C in the initial space is obtained with C=XC1 

\ -0.921916 1.57553 ] 

C = 

-0.0918216 - 1.82313 

with the matrix C, using (4) we build the new matrix P, here it is important to observe the upper 
limit of the summation (n/2), because it implies that the matrix C is trim to the half of columns, 
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P = 2C R C\ , where Cr is a trim matrix of C. This upper limit of the summation implies that n must 
be even, or as is referred in the literature, Restricted Closed Shell. 


P = 2 


-0.921916 

-0.0918216 


[-0.921916 


-0.0918216] 


1.69986 0.169304 

0.169304 0.0168624 


If matrix P is the same as the previous density matrix or within a specified criterion, then 
the procedure has converged. If not, we must calculate again the Fock matrix, using (34), (35) 
and (36) and do all the steps described in this section (MMD). 

Results and Conclusion 

The necessary steps to perform all the matrix process have been programmed in Mathematica 
(M), which is showed in appendix (A), however some particularities of the software must be 
commented. 


• The vectors in M are defined as rows, but the operations with a matrix are with columns. 
Then it is necessary to take the transpose of this vector. 

• M sorts the eigenvalues from highest to lowest, then it is necessary to make additional 
steps to reverse this fact. The solution we found was to build an object with two different 
entities, a scalar as first element and a vector as second element. Then when we calculate 
the transformation matrix ( U ), must be referred to the second element (step 6). But when 
we calculate the unitary matrix ( X ), the denominator is referred to the first element (the 
scalar value in step 8). Again the scalar value was used in step 9. 

• The solution described in the previous item was applied again in step 14, for the 
construction of Fock matrix. 

• In steps 10 and 11 we checked that U and X were correctly budded, with the obtainment 
of S and I. 

• In step 18 we calculate Cr, and this is because the upper limit of the summation is n/2. 

• The Hartree-Fock-Roothaan Energy founded is in agreement with the value reported 
(Szabo, 1996, Levine, 2008). 

• The two-electron repulsion integrals were made with an evaluation formula that we 
construct (17). 

For teaching this concepts in college it is necessary to have enough knowledge in quantum 
mechanics, electromagnetic theory and linear algebra, that is because all the spin orbitals must 
form an orthonormal basis and electron’s repulsion integrals are constructed with those spin 
orbitals. Special attention must have the Hamiltonian that is obtained with the Schrodinger 
equation. 
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Appendix 1: Mathematica Programming 

(*1.- Clear Global memory *) 
ClearAll["Globar* M ] 

(*2.- Define the matrix Ef ore *) 

F1C={ {-1.8488,-1.8826}, {-1.8826,-1.5860}}; 
" H° ore = ”MatrixForm[HC] 


(*3.- Define the overlapping matrix S *) 

S={{1,0.8366}, {0.8366,1}}; 

"Overlapping matrix S = M MatrixForm[S] 

(*4.- Define an object that is composed of a scalar (eigenvalue) in the first column and a vector 
(eigenvector) in the second column, which are obtained by diagonalizing the matrix S *) 
{eigenvals,eigenvecs}=Eigensystem[S]; 

(*5.- Assign to a variable 'so' the object organized from lowest to highest *) 

"Sorted object that was obtained from the diagonalization of the overlapping matrix" 
so=SortBy[Transpose[{eigenvals,eigenvecs}], First] 

(*6.- Take the element in the second row of the organized object 'so' that are the eigenvectors and with 
these form a matrix *) 

MMS={so[[l,2]],so[[2,2]]}; 

"Eigenvectors of the sorted object defined as rows = "MatrixForm[MMS] 

(*7.- The transformation matrix U, is formed with the transpose of MMS (the eigenvectors are written as 
column vectors)*) 

U=Transpose[MMS]; 

"U = "MatrixForm[U] 

(*8.- The transformation matrix X is obtained by dividing the eigenvector's columns of U between its 
corresponding eigenvalue's square root *) 

X={ {U[[ 1, l]]/Sqrt[so[[ 1,1]]],U[[ 1,2]]/Sqrt[so[[2,1]]]}, {TJ[[2, l]]/Sqrt[so[[ 1, l]]],TJ[[2,2]]/Sqrt[so[[2,1]]]}}; 
"Transformation Matrix X = ”MatrixForm[X] 

"Matrix Transpose of X = ”MatrixForm[Transpose[X]] 

(*9.- The diagonal matrix DS is formed with the eigenvalues of S *) 

DS={{so[[1,1]],0},{0,so[[2,1]]}}; 

"Diagonal Matrix of eigenvalues of DS = ”MatrixForm[DS] 

(*10.- Proof that the Matrix multiplication of U with DS and with the transpose of the matrix U returns us 
S*) 

" Proof U.DS.U + = "MatrixForm[U.DS.Transpose[U]] 

(*11.- Proof that the Matrix multiplication of X with DS and with the transpose of the matrix X returns us 
the Unitary matrix *) 

" Proof X.DS.X t = "MatrixForm[X.DS.Transpose[X]] 


(*12.- Define the initial matrix P *) 
p=»U},{U}}; 

"P = "MatrixForm[P] 
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(*13.- The loop to find the lowest energy is initiated *) 

For[i=l,i<=7,i++, 

Frs={{HC[[l,l]]+0.4531P[[l,l]]+0.9033P[[l,2]]+0.7058P[[2,2]],HC[[l,2]]+0.45165P[[l,l]]+0.8391P[[l, 

2]]+0.6490P[[2,2]]},{HC[[2,l]]+0.45165P[[l,l]]+0.8391P[[l,2]]+0.6490P[[2,2]],HC[[2,2]]+0.7058P[[l,l 

]]+1.2980P[[l,2]]+0.9094P[[2,2]]}}; 

Print[ M Fock Matrix F rs ] = " MatrixForm[Frs]]; 

Fl=Transpose[X].Frs.X; 

Print[ ?, Transformed Fock Matrix FI = " MatrixForm[Fl]]; 

(*14.- Diagonalize the Fock matrix transform FI and sort the eigenvalues and their eigenvectors 
associated, from lowest to highest *); 

Clear[eigenvals,eigenvecs] 

{eigenvals,eigenvecs} =Eigensystem[F 1 ]; 
so=SortBy[Transpose[{eigenvals,eigenvecs}], First]; 

(*15.- matrix NF1 is formed with the eigenvectors sorted, taken from object ? so ? with the eigenvectors as 
rows *); 

NFl={so[[l,2]],so[[2,2]]}; 

(*16.- The matrix of coefficients Cl is formed with the transpose of NF1 (the eigenvectors are placed as 
columns)*); 

C1 =T:ranspose [NF1 ]; 

Print[ ?f Matrix of coefficients in orthonormal espace Cl = " MatrixForm[Cl]]; 

(*17.- Place the matrix of coefficients Cl in the initial vectorial space, the matrix C is formed *); 
CM=X.C1; 

Print[ ?f Matrix of coefficients C = " MatrixForm[CM]]; 

(*18.- The elements values of the matrix P are calculated again, with the trimmed matrix C *); 

CR=Transpose [Transpose [CM] [ [ 1;; 1 ] ] ]; 

Print[ ?, Trimmed matrix CR = ”MatrixForm[CR]]; 

P=2CR.Transpose[CR]; 

Print["P = ?, MatrixForm[P]]; 

(*19.- The Hartree-Fock-Roothaan Energy is calculated for the new values of matrix P *); 

EHF 

=0.5(P[[l,l]](Frs[[l,l]]+HC[[l,l]])+P[[l,2]](Frs[[l,2]]+HC[[l,2]])+P[[2,l]](Frs[[2,l]]+HC[[2,l]])+P[[2, 

2]](Frs[[2,2]]+HC[[2,2]])); 

Print["Hartree-Fock-Roothaan Energy = ",EHF]; 

Print[Style[{"Finish Loop ",i},FontColor -> Red]]; 

] 

(*20.- Displays the last energy calculated *) 

Print[ ?f Optimal energy of the system E HF = ", EHF]; 
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